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We study the kinetics of assembly of two plates of varying hydrophobicity, including cases 
where drying occurs and water strongly solvates the plate surfaces. The potential of mean 
force and molecular-scale hydrodynamics are computed from molecular dynamics simula- 
tions in explicit solvent as a function of particle separation. In agreement with our recent 
work on nanospheres [J. Phys. Chem. B 116, 378 (2012)] regions of high friction are found 
to be engendered by large and slow solvent fluctuations. These slow fluctuations can be due 
to either drying or confinement. The mean first passage times for assembly are computed 
by means of molecular dynamics simulations in explicit solvent and by Brownian dynam- 
ics simulations along the reaction path. Brownian dynamics makes use of the potential 
of mean force and hydrodynamic profile that we determined. Surprisingly, we find rea- 
sonable agreement between full scale molecular dynamics and Brownian dynamics, despite 
the role of slow solvent relaxation in the assembly process. We found that molecular scale 
hydrodynamic interactions are essential in describing the kinetics of assembly. 



I. INTRODUCTION 

Hydrophobic interactions are important in 
structural biology and nanoscience where they 
are often the major driving force in self- 
assemblyP^ F r example, hydrophobic inter- 
actions are acknowledged to play a major role 
in the formation and folding of proteins! 10 * 11 ! 
A classic example of hydrophobic assembly is 
when two plates that bear little attraction to 
water are separated by less than some critical 
separation, capillary evaporation occurs in the 
inter-plate region, thereby driving association. 

Although a great deal of work has concen- 
trated on the structural and thermodynamic as- 
pects of hydrophobic assembly, less is known 
about the kinetic mechanisms involved in the 
association process. Prior studies have inves- 
tigated the rate of evaporation of solve nt co n- 
fined between two hydrophobic surfaces EE2KS1 
although the full association pathway remains 



largely unexplored. One possible route to mod- 
eling the kinetics is within the framework of 
Brownian (Smoluchowski) dynamics, a coarse 
grained description in which the solvent degrees 
of freedom are not explicitly treated and where 
the dynamics of the heavy bodies are instead 
modeled in a stochastic bath. The effect of sol- 
vent is encoded in the potential of mean force 
(PMF) between bodies and the friction coeffi- 
cient. If hydrodynamic interactions (HI) are in- 
cluded, then the friction coefficients depend on 
the positions of the heavy bodies. Within this 
framework, one can simulate protein diffusion 
and association using the appropriate equations 
together with the potential of mean force and 
the (position dependent) friction tensor as in- 
put!^ For example, such an approach has been 
utilized to study crowded cellular environments, 
where the diffusion of proteins is thought to slow 
down considerably due in part to hydrodynamic 
interactions.^ 
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Most typically, hydrodynamic interactions 
that are utilized in Brownian dynamics simula- 
tion are computed within a continuum descrip- 



2 



tion of the solvent E§li2l This description, how- 
ever break s dow n for inter-particle length scales 
of 1-2 nmPSEB On such length scales, fluctua- 
tions which accompany solvent layering or cap- 
illary drying could have an enormous effect on 
hydrodynamic interactions. In recent work, we 
have shown how the nature of solvent confined 
between nanoscopic spheres determines the role 
of HI in assembly. 2 2 There is a growing interest 
in the behavior of water molecules in confined 
geometries The diffusion constant and the 
hydrogen bond lifetime of water molecules in 
confined environments are distinct from those 
in the bulkP^l Furthermore, similar signatures 
have been recently observed in simulations of 
crowded protein solutions.^ We believe such 
effects can potentially play an important role 
in the transport and association of nanoscopic 
bodies. 

In this paper, we extend our previous study 
of nanospheres to the study of hydrodynamic 
interactions and association kinetics of plates 
of varying hydrophobicity. We calculate the 
friction on two parallel plates as a function of 
their separation from molecular dynamic sim- 
ulations with explicit water using a technique 
previously employed for spherical solutes 1221281 
We also calculate the potential of mean force 
between the two plates as a function of their 
separation. Both the molecular-scale effects of 
hydrodynamic interactions and the free energy 
profile are found to be highly correlated with the 
behavior of water molecules in the inter-plate 
region, in agreement with our prior work!22l 

We presently consider three types of 
graphene-like plates: plates with a "fully" 
attractive solute-solvent interaction potential, 
plates with a "reduced" attractive solute- 
solvent potential, and plates with a purely 
repulsive solute-solvent potential. For the 
plates with either a reduced attractive poten- 
tial or a purely repulsive potential, capillary 
evaporation occurs when the inter-plate separa- 
tion is smaller than a critical value. In contrast, 
capillary evaporation does not occur for plates 
with a fully attractive solute-solvent potential. 
As the fully attractive plates approach each 
other, we observe solvent layering with water 



molecules being eventually expelled from the 
inter-plate region due to steric repulsion. Both 
dewetting in the former cases and solvent 
layering in the latter case greatly affect the 
solvent fluctuations in the inter-plate region 
and give rise to molecular scale effects in the 
hydrodynamic interactions. We compute the 
spatial dependence of the friction tensor and 
investigate the relation between the static and 
dynamic fluctuations of water density in the 
inter-plate region and this property. 

Given the potential of mean force and posi- 
tion dependent friction coefficient we compute 
the rate of diffusion controlled association by 
means of Brownian dynamics (BD) simulations. 
The predictions of BD are tested against the re- 
sults garnered from a set of molecular dynamics 
(MD) simulations of assembly in explicit sol- 
vent. This serves as a test of the validity of the 
Markovian approximation inherent in Brownian 
dynamics simulation, and the ability of coarse- 
grained stochastic dynamics to adequately cap- 
ture the kinetic effects associated with the bath. 
We find that there is reasonable agreement be- 
tween the two simulations when molecular scale 
hydrodynamic interactions are included in the 
BD, but marked disagreement when hydrody- 
namic interactions are neglected. The observa- 
tion that BD with hydrodynamic interactions 
is in agreement with the MD result is some- 
what surprising since, as we show, slow solvent 
fluctuations play an important role in the as- 
sembly process, and we would thus expect non- 
Markovian effects to be of importance. The 
agreement between BD and MD indicates that 
the non-Markovian effects may be encoded in 
the spatial dependence of the friction coefficient 
in an average way. 



II. SYSTEM AND METHOD 

We modeled the hydrophobic plate as a 
graphene-like sheet with an area of (1.2 x 1.3) 
nm 2 as shown in Fig. 1. The plate size was cho- 
sen so as to facilitate both nanometer scale dry- 
ing transitions and comparison with our prior 
work on nanospheres. In order to vary the hy- 



3 



drophobicity we studied three types of plates, 
all with a carbon-carbon bond length of 0.145 
nm, and a Lennard- Jones diameter of occ = 
0.35 nm but with different interaction potential 
well depths: (a) full Lennard- Jones (LJ) inter- 
action with ecc = 0.276 kJ/molp^ (b) reduced 
LJ interaction with ecc — 0.055 kJ/mol, and 
(c) purely repulsive Weeks- Chandler- Anderson 
(WCA) truncation of the full LJ potential 
For convenience, these three types of plates are 
denoted as (full) LJ plates, reduced LJ plates 
and WCA plates in the present work. Solute- 
solvent interactions are given by the geometric 
mean of the respective water and solute param- 
eters. 

In order to calculate the relative friction co- 
efficient, two parallel plates are placed perpen- 
dicular to the z-axis and a series of simulations 
are run with the plates fixed at various sepa- 
rations, and all internal degrees of freedom are 
frozen. The inter-plate separation ranges from 
0.4 to 2 nm (for the full and reduced LJ plates) 
or to 2.2 nm (for the WCA plates). The solvent- 
induced mean force is computed from such fixed 
configurations, although a finer grid spacing in 
z is employed. In order to obtain the poten- 
tial of mean force (PMF) , this quantity is then 
integrated in combination with the direct plate- 
plate interactions. 

In the simulations of the association process, 
two plates are initially placed f=a 1.8 nm apart. 
The lower plate is fixed to its position and a 
biased force along the — z direction is applied 
to the upper plate. Harmonic potentials are 
utilized to treat the intramolecular stretches 
(k b = 392721.8 kJ/mol nm- 2 and d = 0.14 
nm) and bends (kg = 5271.84 kJ/mol rad -2 
and 8 = 120°). Snapshots from the association 
of reduced LJ plates and representative trajec- 
tories plotted along the reaction coordinate are 
shown in Fig. 1. The results from these two sets 
of calculations are discussed in Section llV Bl and 



Section [IV C[ respectively. In addition, we de- 
termine the friction coefficient on a single plate, 
in a system containing only one plate, by pulling 
it at a constant speed perpendicular to its plane 
and computing the drag force along its normal 
direction. 



All systems are solvated by TIP4P water.^ 
The size of the solvation box of the full LJ and 
reduced L J system is 4 x 4 x 4 nm 3 , and for the 
WCA system the box size is increased to 5 x 5 x 6 
nm 3 . This is done in order to accommodate the 
large volume fluctuation that arising from the 
strong dewetting transition. 

In the calculation of the friction, each system 
is equilibrated with a 2 ns NPT simulation (P= 
1 atm, T = 300 K). The Berendsen method 1 ^ 
and the stochastic velocity rescaling methocP^ 
are chosen for the barostat and thermostat, re- 
spectively. For each system, there are 12-18 
NVE production runs of 4 ns to compute the 
friction coefficient. In order to facilitate energy 
conservation, double precision routines are uti- 
lized for all NVE runs. To study the dragging 
force of the single plate, we choose four pulling 
rates ranging from 0.5 to 3 nm/ns. There are at 
least 12 NPT runs for each pulling rate that are 
utilized to obtain an accurate estimate of drag- 
ging force. We perform at least 55 runs with 
various initial configurations for each system to 
fully characterize the association process of two 
plates and estimate the distribution of mean 
passage times. All simulations are performed 
using GROMACS 4.5.3 and the Particle-Mesh 
Ewald technique is utilized to treat long range 
electrostatics ! 34135 1 



III. BROWNIAN DYNAMICS WITH 
HYDRODYNAMIC INTERACTIONS 

Within the framework of Brownian dynam- 
ics, hydrodynamic interactions are encoded in 
the frictional force. The stochastic equation for 
the association process contains contributions 
from the mean force, the frictional force, and 
the Gaussian random force. The potential of 
mean force, W(z), includes the contributions 
from the direct plate-plate interaction and the 
solvent-induced interaction. The separation be- 
tween two plates, z, is treated as the reaction co- 
ordinate. Hydrodynamic interactions give rise 
to the spatial dependence of the relative friction 
coefficient, £(z) and through this the spatially 
dependent diffusion coefficient D(z) = kT/((z). 
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FIG. 1. (Upper left panel) The top view and side view of two parallel graphene-like plates. (Upper 
right panel) Representative association trajectories of two plates with the "reduced" LJ interaction. (Lower 
panels) Configurations taken from the association of the reduced LJ plates, corresponding to an inter- 
plate separation of (from left to right): z = 1.7 nm (initial configuration), 1.0 nm (where the association 
stagnates), 0.7 nm (during the driving induced collapse) and 0.35 nm (when the plates come into contact). 
The plates and water molecules in the inter-plate region are depicted by a space filling representation, and 
the snapshots are rendered with VMDp^ 



The stochastic BD equation forplate associ- 
ation can be therefore written asP3 



-PD(z)^W(z) 



d_ 

dz 



D(z) + R(z,t), (1) 



where f3 = 1/kgT and the diffusion coefficient 
is related to the "random force" R(z, t) through 
the fluctuation dissipation theorem, 



{R(t')R{t)) = 2D(z)5{t-t'). 



(2) 



This equation of motion is an accurate descrip- 
tion of the dynamics in the high friction limit 
for timescales greater than the momentum re- 
laxation time, and in which the solvent time- 
scale is fast compared to the time scale for the 
motion of the heavy bodies. The dynamics can 
alternatively be expressed as a Smoluchowski 
equation for the probability distribution P(z, t) 
of finding the particle at z at time t: 

dP(z,t) d , s ( d A , , 

(3) 

This equation or BD can be used to determine 
mean first passage times for diffusion controlled 
reactions once W(z) and are known. 



In order to apply BD we must determine 
the PMF and spatial dependence of the fric- 
tion coefficient along the reaction path. As 
already pointed out, the continuum treatment 
of the hydrodynamic interaction breaks down 
at small separations pHHUl j n ^ n j s wor k, we cal- 
culate the friction coefficient on the molecular 
length scale from explicit molecular dynamics 
simulations. The friction coefficient can then 
be determined from the two-body friction ten- 
sor Qj. In the calculation, the plates are fixed 
as required in the Brownian limit. The pair fric- 
tion tensor can be expressed as a time integral 
of the correlation functions of the fluctuating 
force 5F i =F i -{F i ), 

C« = P / dt lim (^(^(O)), (4) 

where by symmetry £12 = C21 an d £11 = £22- 
Here we focus on the direction parallel to the 
inter-plate separation and therefore study the 
force components along the plate normal vec- 
tor. Eqn. [4] is only valid in the limit where the 
number of solvent molecules, n, approaches in- 
finity!^! i n the present case of finite systems, it 
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is still possible to relate the force autocorrela- 
tion function to the friction tensor, following the 
procedure of Bocquet et alP^The friction along 
the relative coordinate is given by (£n — Ci2)/2 



and depicted in Section IV B 



IV. RESULTS 

A. Friction and diffusion coefficients of a single 
plate 
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The friction coefficient perpendicular to the 
face of a single plate immersed in solvent may be 
extracted from the force autocorrelation func- 
tion by means of techniques outlined in Ref. 
[571 The resultant calculation yields values of 
1.63X10" 11 , 1.33X10" 11 , l.OxlCT 11 kg/s, when 
the solvent-solute interaction is the full LJ, Re- 
duced LJ, and the WCA potential, respectively. 
In agreement with prior workp^l the friction ex- 
perienced by the body decreases with increasing 
surface hydrophobicity. Each calculation is per- 
formed in a periodic box of dimensions given in 
Section HB 

The friction coefficient may also be extracted 
from a set of non-equilibrium simulations where 
the plate is dragged through the solvent at a 
given velocity. As depicted in Fig. 2, the drag- 
ging force at four velocities is determined from 
MD simulations in explicit solvent, and grows 
with increasing velocity. The resultant plot of 
dragging force versus pulling rate exhibits a lin- 
ear relationship with the slope equal to the fric- 
tion coefficient. The extracted values of the fric- 
tion are given in Table 1, alongside the results 
computed from the force autocorrelation func- 
tion. The two estimates of the friction are in 
close agreement. The consistency between the 
results garnered from two different techniques 
serves to validate the calculation. 

Next, we compare these results to the predic- 
tions of continuum hydrodynamics. The fric- 
tion coefficient for a rigid body described by a 
set of spherical sites can be evaluated by means 
of the Rotne-Prager tensor.^ The translational 
friction coefficient may be extracted from the 
full 37V x 37V mobility tensor, where N is the 



FIG. 2. The dragging force of single plate mov- 
ing at constant velocity. The results for the plate 
with full Lennard-Jones(LJ) interaction (LJ plate), 
the plate with reduced LJ interaction (reduced LJ 
plate), purely repulsive Weeks-Chandler- Anderson 
truncation of the full LJ potential (WCA plate) are 
depicted in black, red, and green, respectively. The 
dragging forces increase linearly as the pulling rate 
increases, and the fitted slopes are equal to the fric- 
tion coefficient. 



number of sitesP^In the present calculation, pe- 
riodic boundary conditions must be taken into 
account and are known to have a significant im- 
pact on hydrodynamic properties! " I Periodic 
effects may be treated by means of replacing the 
standard Rotne-Prager tensor with a form that 
accounts for periodicity via an Ewald summa- 
tionPH This approach is discussed further in the 
Appendix. 

As discussed above, different solvation boxes 
are utilized in the WCA and Lennard- Jones sys- 
tems, and the continuum calculation must be 
undertaken for both sizes. The values of the 
friction coefficient with stick boundary condi- 
tions for the two box sizes are given in Table [ij 
The full LJ value is somewhat larger than the 
continuum result whereas, in line with expecta- 
tions, the increasingly hydrophobic plates begin 
to show larger deviations from the stick bound- 
ary condition. 
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interaction 


box (mil 3 ) 


Cfacp ( 1CT 11 kg/s) 


Cpull ( lO" 11 kg/s) 


CSTICK (HT 11 kg/s) 


WCA 


5x5x6 


1.00 


1.03 


1.27 


REDUCED 


4x4x4 


1.33 


1.36 


1.53 


FULL 


4x4x4 


1.63 


1.69 


1.53 



TABLE I. The friction on a single plate in the z-direction, as estimated from the force autocorrelation 
function, pulling simulations, and continuum hydrodynamics. 



B. Hydrodynamic and thermodynamic profiles 

Friction coefficients, C( z ): an d the potential 
of mean force, W(z), are depicted as a function 
of inter-plate separation, z, in Figs. 3, 4, and 5 
(for the full LJ, reduced LJ, and WCA plate, re- 
spectively). The frictional profiles exhibit non- 
monotonic behavior as the two plates approach 
each other. The spatially dependent features of 
the molecular scale hydrodynamic interactions 
display similar trends to those found in the free 
energy profile. This finding is in agreement with 
our previous work and the recent work of Mit- 
tal and Hummerp21 In the case of the full LJ 
plate, the friction peaks and the free energy in- 
creases as a layer of water is "squeezed out"^ 
at z = 0.88 nm. The friction coefficient subse- 
quently decreases at the minimum of the PMF. 
The final solvent layer is expelled as the sepa- 
ration decreases to 0.62 nm, and both the free 
energy and the friction coefficient increase. It is 
important to note that the friction profile peaks 
at 0.88 and 0.62 nm which are at the same posi- 
tions as the PMF barriers. The peak heights at 
these two separations are 23.3(4.2) x 10~ 10 and 
4.0(0.4) x 10~ 10 kg/s, respectively. Standard er- 
rors of the mean value are given in parenthesis. 

For the reduced LJ plates, there is a low bar- 
rier at z = 0.9 nm in the PMF, whereas the 
WCA plates exhibit barrierless assembly along 
the chosen reaction coordinate. In these cases, 
the corresponding friction profiles also display 
non-monotonic behavior. The friction profile 
peaks at z c = 0.92 nm and z c = 1.35 nm for 
the reduced LJ and WCA plates, respectively. 
As discussed below, these distances are in the 
region of the critical separation for dewetting. 
The corresponding peak heights of the friction 
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FIG. 3. Spatial dependence of relative friction co- 
efficient (A) and the potential of mean force (B) for 
the association of two L J plates. The friction peaks 
at z = 0.62 nm and two free energy wells at z — 
0.68 and 0.97 nm are depicted in the insets of (A) 
and (B), respectively. (C) The relative density of 
water in the inter-plate region. (D) The ratio of 
the variance to the average of the number of water 
molecules (black) and the solvent relaxation time in 
the inter-plate region (red). For both curves in (D), 
the peaks at z « 0.9 nm are depicted in the inset. 
The values corresponding to the friction peaks are 
shown as filled symbols in panels (C) and (D). 



coefficients are 4.1(0.85) x 10" 10 (reduced LJ) 
and 2.2(0.61) x lO" 10 kg/s (WCA). The some- 
what large error bar results from the sizable 
force fluctuations that are present in the vicin- 
ity of z c . 

The friction profiles converge to 2.1 x 10 -11 
(LJ), 2.2 x 10" u (reduced LJ) and 1.6 x 10" 11 
(WCA) kg/s at large separations. The standard 
error in the mean for these values is ~ 10%. 
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These numbers are somewhat larger than ex- 
pected upon comparison with the results given 
in Table 1. This difference may largely result 
from the finite box size and the impact of neigh- 
boring periodic images. However, the effect of 
periodic boundary conditions should not greatly 
impact the observed spatial dependence of the 
short-range hydrodynamics interactions. 

In order to better understand the nature 
of molecular-scale hydrodynamic interactions, 
we analyze the density and the static and dy- 
namic fluctuations in number of water molecules 
in the inter-plate region. The static fluctua- 
tions are measured by the ratio of the num- 
ber variance to the average ((SN) 2 )/ (N). An 
estimate of the solvent relaxation time can be 
computed from the integral of the normalized 
autocorrelation function of these fluctuations, 
(SN{t)SN(0))/((SN) 2 ). The resultant plots of 
these quantities as a function of inter-plate sep- 
aration are depicted in Figs. 3-5. Addition- 
ally, the values of the solvent density and fluc- 
tuations at the separation that corresponds to 
the maximum friction are marked on the curves. 
For the LJ plates, the solvent density decreases 
when the plate separation decreases from 0.9 to 
0.85 nm. The ratio of the variance to the aver- 
age of water density peaks at z = 0.88 nm, and 
the relaxation time peaks at z = 0.90 nm. Both 
peaks occur at separations close to where the 
friction coefficient peaks. Moreover, when the 
separation decreases from 0.65 to 0.62 nm the 
solvent density sharply decreases. The static 
fluctuations and relaxation time grow in the 
process of water expulsion. Both the static vari- 
ance of water density and the relaxation time 
peak at the same separation where the friction 
coefficient peaks (z c = 0.62 nm). 

In both the reduced LJ and WCA systems, 
the solvent density in the inter-plate region 
begins to decrease at the critical separation 
of the dewetting transition. For the reduced 
LJ plate, the solvent density dramatically de- 
creases when the separation decreases from 0.95 
to 0.92 nm. Both the variance of water den- 
sity and the solvent relaxation time peak at 
z c = 0.92 nm. The density of water between 
the WCA plates decreases as the separation de- 




z (nm) z (nm) 



FIG. 4. Spatial dependence of relative friction co- 
efficient (A) and the potential of mean force (B) for 
the association of two reduced LJ plates. (C) The 
relative density of water in the inter-plate region. 
(D) The ratio of the variance to the average of the 
number of water molecules (black) and the solvent 
relaxation time in the inter-plate region (red). The 
value corresponding to the friction peak is shown as 
filled symbols in panels (C) and (D). 



creases from 1.45 to 1.35 nm. The ratio of the 
variance to the average of water density and 
the solvent relaxation time peak in the same 
region. At the critical separation for the dewet- 
ting transition, the inter-plate region fluctuates 
between wet and dry. The value of the crit- 
ical distance between surfaces is on the order 
of one nanometer for the present-sized plates 
and decreases with decreasing hydrophobicity, 
in agree ment w ith macroscopic thermodynamic 
analysis! 6 * 15 * 44 ! This characterization of the den- 
sity and the fluctuation of solvent is co nsistent 
with the results of previous studiesP^H 

In the three systems presently studied, it can 
be clearly seen that the friction coefficient in- 
creases where the solvent fluctuations become 
large and slow. Taken together, both static and 
dynamic solvent behavior engender the large 
frictions at small inter-plate separations. In 
agreement with our prior work, 22 molecular- 
scale hydrodynamic interactions largely result 
from such fluctuations when, in the case of hy- 
drophobic bodies, the drying transition occurs 
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FIG. 5. Spatial dependence of relative friction co- 
efficient (A) and the potential of mean force (B) 
for the association of two WCA plates. (C) The 
relative density of water in the inter-plate region. 
(D) The ratio of the variance to the average of the 
number of water molecules (black) and the solvent 
relaxation time in the inter-plate region (red). The 
value corresponding to the friction peak is shown as 
filled symbols in panels (C) and (D). 



or, for less hydrophobic species, when water 
molecules are expelled from the inter-plate re- 
gion due to steric repulsion. 

C. Comparison of Brownian dynamics with 
molecular dynamics 

In order to further elucidate the impact of 
molecular-scale hydrodynamics on the kinetics 
of self assembly, we perform direct molecular 
dynamics simulations of this process. The two 
plates are initially placed perpendicular to the 
z-direction at a separation of ~1.8 nm. A con- 
stant loading force is added to the upper plate 
and the lower plate is fixed in its initial posi- 
tion. The loading force utilized is 440 pN (full 
LJ) or 20 pN (reduced LJ and WCA). There 
are at least fifty-five association simulations per- 
formed for each system. 

The assembly of both the WCA plates and 
reduced LJ plates slows down around the crit- 
ical separation for dewetting. This process is 



illustrated for the reduced LJ ca se by the sam- 
ple trajectories depicted in Fig. |llf^ One can 
see there is a large dwell time at z ~ 1 nm. In 
the case of the full LJ system, the upper plate 
initially rapidly approaches the lower plate and 
the inter-plate distance decreases to z = 0.95 
nm, corresponding to two inter-plate water lay- 
ers. The association process then stagnates. Af- 
ter some period a water layer is expelled and 
the plate separation decreases to z = 0.66 nm 
where a single water layer separates the plates. 
There is another dwell time before the last water 
layer is expelled and the two plates finally come 
into contact. During the association process the 
upper plate can rock. This is particularly pro- 
nounced around z = 0.66 nm where the ampli- 
tude of the rocking can be as much as 0.2 nm. In 
order to provide the best comparison with the 
(one-dimensional) BD scheme outlined below, 
the plates must be kept as parallel as possible. 
This is facilitated by the addition of a harmonic 
restraining potential on the plate's internal de- 
grees of freedom. It is important to note that 
outside the present context, "rocking" could be 
a viable degree of freedom, and we observe that 
assembly occurs significantly more rapidly when 
such effects are included. In the case of less at- 
tractive or purely repulsive solvent-solute inter- 
actions, the effect is much less prominent, as the 
displacement generated by the rocking mode is 
small with respect to the critical separation for 
dewetting that drives assembly (ss 1 nm). 

As we will show, the mean first passage times 
observed in the MD association trajectories can- 
not be predicted by solely considering the free 
energy profile with constant friction. In order to 
evaluate the role of hydrodynamic interactions 
in assembly, we utilize a one-dimensional Brow- 
nian dynamics (BD) scheme as described in Sec- 
tion |III| where the system evolves according to 
Eqn. [T] and can be inte gra ted as described by 
Ermak and McCammoni™ The dynamical de- 
gree of freedom is taken to be the inter-plate 
distance, z, along which the friction and poten- 
tial of mean force have been computed (see Figs. 
3-5). 

This scheme may be utilized to generate pre- 
dictions for the mean first passage time (mfpt) 
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type 


F (pN) 


tmd (ps) 


tbd-hi (ps) 


tbd-nohi (ps) 


z (nm) 


z A (nm) 


WCA 


20 


515 


516 


134 


1.60 


0.40 


reduced LJ 


20 


1469 


2390 


200 


1.20 


0.40 


LJ (1st barrier) 


440 


573 


1460 


85 


1.20 


0.80 


LJ (2nd barrier) 


440 


32600 


23078 


243 


0.75 


0.40 



TABLE II. Mean first passage time, r, of the plate association process as described in Section [IV C| 



and distribution of first passage times (dfpt) 
from BD simulations with the loading force. 
Such a comparison serves to evaluate the degree 
to which the Brownian framework can yield an 
adequate description of the kinetics of assembly. 
The constant force Fq applied to the molecular 
dynamics simulations is accounted in the BD by 
a means of the modified potential of mean force 
U(z), 



U{z) = W(z)+F z, 



(5) 



where W{z) is the PMF determined from MD 
in the absence of a loading force. Because the 
loading force Fq in Eqn. [5] is independent of 
the solution degrees of freedom, Eqn. [5] is rig- 
orous. Moreover, because the loading potential 
is linear in z we assume it will not alter the 
spatial dependence of the friction. However, for 
very large loading forces we expect that non- 
equilibrium effects will be significant and the 
BD picture will fail. For each system, we have 
generated 10,000 BD trajectories given initial 
state zq, and an absorbing boundary at za- For 
each system, the values of these parameters and 
of the mean first passage time are given in Table 
[TT| The corresponding first passage time distri- 
butions are shown in Fig. 6. The mean first pas- 
sage time may also be computed directly from 
the Smoluchowski equation (Eqn. [3]) by means 
of the following expression!^ 



T ( z o) = P dy 



C(y) 

-0U( V ) 



dxe-P u< - x \ (6) 



2R 



where f3 = l/fc^T and zr is the positon of 
the reflecting boundary. Due to the the driving 
force, the reflecting boundary can be taken to 



large values in our calculation. This equation 
and an ensemble of BD trajectories will yield 
the same result within numerical error. 

The mean first passage time obtained from 
MD simulation for the association of the re- 
duced LJ plates, is about 60% smaller than the 
result from BD simulation with hydrodynamic 
interaction (BD with HI), but 7 times larger 
than the value estimated from BD when the spa- 
tial dependence of the friction is not included 
and the single plate friction (see Table |l| is uti- 
lized instead (BD without HI). The fpt distri- 
bution obtained from MD is also close to, albeit 
narrower than, that obtained from BD with HI. 
In contrast, the distribution garnered from BD 
without HI is much more strongly peaked. In 
the case of the WCA plates, the mfpt and the 
fpt distribution from both MD and BD with HI 
are very close to each other, but both are ap- 
proximately 4 times larger than the BD without 
HI result. 

The results of MD simulation are in reason- 
able agreement with those obtained from BD 
with HI for both the WCA and the reduced LJ 
plate systems. The association process slows 
down in the region around the critical separa- 
tion for dewetting transition where the friction 
peaks. As discussed above, the behavior of fric- 
tion profile at the critical separation largely re- 
sults from the solvent fluctuation due to the 
dewetting transition. The molecular-scale ef- 
fect of hydrodynamic interaction evidently con- 
tributes to the slowing down of the association 
process near the critical separation, and if only 
barriers present in the PMF are considered (as 
in BD without HI), association occurs far too 
rapidly. These results indicate that a kinetic 
barrier along the reaction coordinate is present 
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FIG. 6. The distribution of first passage time (dfpt) 
for the association process of LJ plates (panel (A), 
first barrier; panel (B), second barrier), reduced LJ 
plates (C) and WCA plates (D), obtained from the 
MD simulation (black), BD simulation with and 
without the consideration of hydrodynamic inter- 
actions (HI) (red and blue, respectively). The df- 
pts obtained from MD are in reasonable agreement 
with the result of BD with HI, and the probability 
of having a fpt smaller than 500 ps is much larger 
in the results garnered BD without HI. 



at the drying transition. 

The values of the mfpt of the first (second) 
barrier of the full LJ plate system as obtained 
from BD with HI are 1460 (23078) ps, much 
larger than the results of 85 (243) ps for BD 
without HI. Hydrodynamic interactions strik- 
ingly slow down the first passage time over the 
second barrier by about two orders of magni- 
tude. The mfpt obtained from the MD sim- 
ulation is 573 and 32600 ps, for the first and 
second barrier, respectively. For the first bar- 
rier, the mfpt obtained from MD simulations is 
smaller than the result from BD with HI, while 
still being much larger than the result from BD 
without HI. Meanwhile, for the second barrier, 
the mean first passage time from MD simula- 
tions is larger than the result from BD sim- 
ulations with HI. The distribution of the first 
passage times corresponding to the second bar- 
rier obtained from MD simulations is similar to 



the results from BD simulations with HI. For 
passage over the first barrier, the distributions 
exhibit greater deviation. 

In general, the average value and the distri- 
bution of mean first passage times of the three 
types of plates calculated from MD simulations 
is consistent with the results of BD simula- 
tions with HI. In prior workp^l we found that 
hydrodynamic interactions contributed approx- 
imately 40% to the assembly of two fullerenes. 
In the present set of calculations, hydrodynamic 
interactions contribute a much larger share. 
Comparison of spheres and plates indicate that 
the contribution of the hydrodynamic interac- 
tion is enhanced as the shape is flattened. More 
water molecules are confined by plates than 
spheres of the same surface area so that, at 
small separations, the degree of confinement 
and the length scale of dewetting is increased. 



V. CONCLUSION 

In this work we study the impact of hydrody- 
namic effects on the kinetics of assembly of two 
plates of varying hydrophobicity. To this end, 
the potential of mean force and spatially depen- 
dent friction coefficient are determined along 
the inter-plate separation. The results show 
that there is a correspondence between peaks 
in the PMF and peaks in the frictional profile. 
High values of the friction are related to large 
and slow solvent fluctuations in the inter-plate 
region. Both solvent confinement and drying 
phenomena can play a critical role in the kinet- 
ics of assembly. In this way, molecular scale ef- 
fects shape the hydrodynamic interactions, and 
give rise to their deviation from continuum the- 
ory, which predicts that the friction coef ficien t 
diverges as two bodies come into contact HZBH 

The kinetics of assembly studied by means 
of molecular dynamics simulations in explicit 
solvent with a constant loading force applied 
along the reaction coordinate is compared with 
the predictions of Brownian dynamics with the 
same loading force and with the potential of 
mean force and hydrodynamic profile extracted 
from the data presented in Section |IVB| Brow- 
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man dynamics is a widely used technique ow- 
ing to its computational efficiency as solvent 
degrees of freedom are not treated explicitly. 
There is reasonable agreement between the 
mean first passage times that are obtained from 
the two schemes. Indeed, we find that hydro- 
dynamic interactions are essential to produce 
a reasonable description of the process and ne- 
glect of spatial dependence of the friction has a 
large impact on the kinetics. The HI give rise 
kinetic bottlenecks along the association path- 
way,^ over and above the barriers in the po- 
tential of mean force. Interestingly, other de- 
grees of freedom such as plate "rocking" can 
in some cases significantly increase the rate of 
assembly, probably by facilitating waters entry 
and exit from the inter-plate gap. In order to 
describe the effect of solvent relaxation on the 
rocking mode within the Brownian framework, 
it would require a determination of the friction 
coefficient experienced on this degree of free- 
dom. The exploration of all possible pathways 
to plate assembly is beyond the scope of this 
work, and presently we concentrate on the ap- 
proach of two parallel plates. 

It has been suggested^ that including sol- 
vent degrees of freedom in the reaction coor- 
dinate is necessary in order to properly char- 
acterize the association process. Presently, a 
reaction coordinate is utilized that is only a 
function of the plate degrees of freedom. We 
have shown that such a choice is reasonable 
provided that (molecular-scale) hydrodynamic 
interactions are considered along the pathway. 
The inclusion of hydrodynamic effects captures 
the effect of solvent in an average sense, and the 
Brownian framework would be an exact treat- 
ment in the limit where the solvent time scales 
are much faster than those associated with the 
heavy bodies. Leading from the last point, the 
discrepancies between the MD and BD results 
can be largely attributed to cither the impact of 
the pulling force, that is the friction experienced 
along the force-biased surface significantly dif- 
fers from that extracted from the equilibrium 
calculations, or to non-Markovian effects. In 
the case of the passage over the first barrier 
of the LJ system, at least some of the devia- 



tion is likely due to the large pulling force, al- 
though non-Markovian effects should have an 
impact both in this case and in the other sys- 
tem studied. It has, for instance, been shown 
that dynamical caging effects can be exhibited 
along degrees of freedom that are associated 
with slowly varying memory functions.^ Such 
phenomena can also serve to explain the long 
dwell times exhibited in the association trajec- 
tories. 

Unfortunately, given the broad distributions 
involved and limited number of MD trajecto- 
ries that can be harvested, it is difficult to 
fully gauge the impact of non-Markovian effects 
based on the present work. Certainly, it is ob- 
served that, in the case of strongly hydrophobic 
plates, assembly is almost always preceded by 
the "first" drying transition, that is the sys- 
tem does not sample successive wet and dry 
states near the critical separation. One should 
expect such a process to be intrinsically non- 
Markovian, although the BD with HI approach 
that is presently employed appears to at least 
capture some aspect of this "waiting period" for 
drying in an average way. A more precise under- 
standing of this phenomena will be the subject 
of future work. 
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Appendix A: Continuum hydrodynamic treatment 
of the friction on a rigid body 

Consider a rigid body made up of N spherical 
sites where the origin is taken to be the center 
of mass. A 3N x 3iV mobility tensor B may be 
defined where: 

V = B«F, (Al) 
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and where V and F are 3-/V dimensional vectors 
representing the velocities and external forces 
on the N sites that comprise the body. In the 
case of stick boundary conditions, the elements 
of B are approximated by the Rotne-Prager ten- 
sor^ where the 3x3 elements 



as: 



B . . , are given 



where a is the site radius, and rj is the shear vis- 
cosity. For TIP4P water the value is 77 = 0.494 
mPa s and has been taken from the literature!^] 
The matrix T. . is given by the following expres- 



B. = 7 

— l J onr/a 



£y + (l-%)T. 



(A2) 




where = (ryl and ff.. is the vector direct 
product of the unit vector of displacement. If 
the plate is treated as a body that is centro- 
and axi- symmetric, then the translational fric- 
tion tensor, £ , is given by the following expres- 
sion!^ 

y 

where C is a diagonal 3x3 matrix. Presently 

we are interested in the friction coefficient asso- 
ciated with motion perpendicular to the face of 
the plates, and this is what is reported in Sec- 
tion [TV^] For periodic systems, the elements of 
B are given by 

B pbc = y- B a (A5) 

n 

where n is the periodic image vector, and 
the same as in Eqn. |A2| except it is evaluated 
over periodic images. This expression may be 
evaluated by means of an Ewald summation, as 
was shown by BeenakkerP^ In the case of over- 
lapping spheres, reciprocal space contributions 
are excluded from the Ewald sum, as the long- 
range portion of the hydrodynamic interaction 
only includes terms for which > 2a. 
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